Convective Heat Transport in Compressible Fluids 
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We present hydrodynamic equations of compressible fluids in gravity as a generalization of those 
in the Boussinesq approximation used for nearly incompressible fluids. They account for adiabatic 
processes taking place throughout the cell (the piston effect) and those taking place within plumes 
(the adiabatic temperature gradient effect). Performing two-dimensional numerical analysis, we 
reveal some unique features of plume generation and convection in transient and steady states of 
compressible fluids. As the critical point is approached, overall temperature changes induced by 
plume arrivals at the boundary walls are amplified, giving rise to overshoot behavior in transient 
states and significant noises of the temperature in steady states. The velocity field is suggested to 
assume a logarithmic profile within boundary layers. Random reversal of macroscopic shear flow 
is examined in a cell with unit aspect ratio. We also present a simple scaling theory for moderate 
Rayleigh numbers. 

PACS numbers: 44.25. +f, 47.27.Te, 64.70.Fx 



I. INTRODUCTION 

Recently much attention has been paid to organized 
fluid motion in turbulent convection in the Rayleigh- 
Benard geometry Although the conventional hy- 

drodynamic equations are constructed for (nearly) in- 
compressible fluids H , we may mention a number of con- 
vection experiments in compressible one-component flu- 
ids in the supercritical region [|g| — p~9|] , together with those 
in non-critical fluids such as water or Hg |2(]-[2^]. In 
these studies the Nusselt number Nu representing the ef- 
ficiency of convective heat transport has been measured 
at large values of the Rayleigh number Ra defined by 



Ra = a p pgL 3 AT/r)D. 



(1.1) 



Here g is the gravity constant, AT = Tbot — Ttop is the 
difference between the bottom and top temperatures, 
and L is the cell height. As the critical point is ap- 
proached in one-component fluids, the thermal expansion 
coefficient a p — — (dp/8T) p /p grows strongly as £?l v (in 
the same manner as the isothermal expansion coefficient 
Kt — (dp/dp)x/p and the isobaric specific heat C p ), 
the thermal diffusivity D decreases as and the shear 
viscosity rj is nearly a constant. Here £ is the thermal 
correlation length growing as (T/T c — on the crit- 
ical isochore with 7 = 1.24 and v = 0.625. Hence, in 
the critical region, Ra can be extremely large; for exam- 
ple, Ra ~ 10 13 even for not very long L(< 10 cm). The 
Prandtl number Pr = rj/ pD was in the range of 1-100. 

High compressibility of supercritical fluids gives rise 
to some unique features not encountered in incompress- 
ible fluids, (i) First, the transient behavior after appli- 
cation of a heat flux from the bottom is strongly in- 
fluenced by the so-called piston effect p^-|3C[|, as re- 
vealed by recent high-precision experiments on 3 He 
and reproduced by subsequent simulation (ii) Sec- 

ond, as a p grows, the usual mechanism of convection 



onset Ra > Ra c (= 1708) is replaced by that of the 
Schwarzschild criterion |)^,[33|. That is, for large com- 
pressible fluid columns (even far from the critical point), 
convection sets in when thermal plumes continue to rise 
upward adiabatically. This occurs when the applied tem- 
perature gradient \dT/dz\ is larger than the adiabatic 
gradient Q, 



a g = (dT/dp) s pg, 



(1.2) 



which is equal to 0.034 mK/cm for 3 He and 0.27 mK/cm 
for CO2. This is the condition that the entropy s = 
s(T,p) per unit mass decreases with height as ds/dz = 
(C p /T)[dT/dz + a g ] < 0, under which the entropy of 
fluid elements adiabatically convected upward is larger 
than that of the ambient fluid. More precisely, Gitter- 
man and Steinberg J3^] found that the convection onset 
for compressible fluids is given by Ra colr > Ra c , where 
Ra colT is a corrected Rayleigh number defined by 



Ra c 



(a p pgL 3 /r)D)(AT-a s L) 
Ra(l-a s L/AT). 



(1.3) 



This is a natural consequence because the effective tem- 
perature gradient seen by the raising plumes is given by 
AT/L — cig. At the convection onset we thus have 



(AT), 



Ra c Drj / (gpa p L 3 



(1.4) 



where the second term behaves as (T/T c — \) 1+V L 3 
and can be much smaller than the first term even for 
small L(~ 1mm Q) as T -> T c . The relation (1.4) 
has been confirmed in SF 6 pi, and in 3 He p!q |. (iii) 
Third, in steady convective states, experimental curves 
of Ra(Nu — 1) vs Ra co " were collapsed onto a single 
universal curve for various densities above T c Jl2[ and 
for various average reduced temperatures on the critical 



1 



isochore |L6| . These empirical results are highly nontriv- 
ial, because Nu can in principle depend on Ra, Pr, and 
a g L/AT, while Nu is a function of Ra and Pr for in- 
compressible fluids. 

For various fluids under relatively large AT 3> a g L 
(where Ra corr = Ra), data of Nu have been fitted to a 
simple scaling law, 



Nu ~ Ra a . 



(1.5) 



The exponent a has been in a range from 0.28 to 0.31 
and, in particular, a theoretical value 2/7 |2,p| was gener- 
ally consistent with data for Ra < 10 12 |"l8|,|l). More- 
over, measurements of the patterns of isothermal surfaces 
po| and the velocity [^l|,^2| have been informative on 
plume motion and a large-scale circulating shear flow in 
small- aspect-ratio cells |l^,^l[. Several authors have also 
performed numerical analysis of convection at large Ra 
in two dimensions (2D) p5|-^9| and in three dimensions 
(3D) p?j|-^3|. Even in 2D salient features in the experi- 
ments have been reproduced. In these simulations, if the 
temperature is averaged over a long time, the temper- 
ature gradient is localized in thin boundary layers with 
thickness It related to Nu by 



Nu = L/2£ 7 



(1.6) 



Both in 2D and in 3D (if visualized from side), the plumes 
tend to be connected from bottom to top for large Pr be- 
cause of slow thermal diffusion, while they become diffuse 
far from the boundaries for small Pr. In the 3D simu- 
lations with periodic or free-slip sidewalk |4^,^3|], local 
boundary shear flows were observed between incoming 
plumes and outgoing networks of buoyant sheets in hor- 
izontal planes close to the boundaries. 

In this paper we will derive and examine hydrodynamic 
equations for compressible fluids under gravity in the 
supercritical region, in which the oscillatory motion of 
sound has been averaged out (23) . Since the time scale of 
convective motions is much longer than that of the acous- 
tic wave i ac = L/c (typically of order 10 _4 s for L ~ 1 
cm), such a description is convenient theoretically and is 
even indispensable for numerical analysis. Our dynamic 
equations are a natural generalization of the usual hydro- 
dynamic equations [||. Our new predictions are unique 
particularly when the piston effect comes into play, as 
has been demonstrated in the previous simulation |3l| l 
relatively close to the convection onset. This paper will 
present 2D simulations of our hydrodynamic equations 
for much larger Ra co " both in transient and (dynamical) 
steady states. Even in steady states, we will find some 
characteristic features of turbulent states, which have not 
been reported in the previous simulations [|35[- ^7| , ^0[ -^3| , 
such as the logarithmic velocity profile of the velocity 
near the boundary p4j and random reversal of the large- 
scale circulating flow in small-aspect-ratio cells 
We will also point out that individual arrivals of plumes 



at the boundaries cause global temperature fluctuations 
in the cell via the piston effect. The resultant noise level 
of the temperature fluctuations grows as the critical point 
is approached. 



II. THEORETICAL BACKGROUND 

A. Hydrodynamic equations 

We consider a supercritical fluid on the critical isochore 
in a cell with the bottom plate at z = and the top plate 
at z = L. The z axis is taken in the upward direction 
and the total fluid volume is fixed at V. The temperature 
disturbance ST(r, t) = T(r, t) — T top measured from the 
temperature T top at the top boundary is much smaller 
in magnitude than T top — T c . Hereafter e will be used 
to denote the reduced temperature at the top boundary, 
which satisfies 



e = T top /T c - 1 » AT/T C 



(2.1) 



We assume that the gravity-induced density stratification 
is not too severe such that the thermodynamic derivatives 
are nearly homogeneous in the cell. This is satisfied when 
\p/p c - 1| ~ {dp/dp) T gL < e 13 with 13 = 0.33 @. This 
condition is rewritten as 



J3+7 



> a g L/T c 



(2.2) 



In the theoretical literature on convection |^| the 
top and bottom temperatures Tbot and T top are constant 
parameters. However, in most of recent convection ex- 
periments, particularly in cryogenic ones, the heat flux 
at the bottom Q = —\(dT/dz) z= o and T top have been 
fixed. The A is the thermal conductivity. Furthermore, 
if the top and bottom walls are made of a metal with 
high thermal conductivity, the boundary temperatures 
become homogeneous in the lateral directions (unless lo- 
cal temperature changes are too fast). Then Tbot(i) and 
hence AT(i) are functions of time only. Metcalfe and 
Bchringer [^5| performed linear stability analysis of con- 
vection onset under this cryogenic boundary condition. 
In the nonlinear regime it is of great interest how the 
boundary layer thickness and the plume generation de- 
pend on the boundary condition. 

In equilibrium the pressure gradient is given by — pg = 
—p c g. In noncquilibrium we set 



p(r,t) =p - p c gz + px{t) +Pinh(r,t), 



(2.3) 



where po is a constant, pi(t) and pi n h are the homoge- 
neous and inhomogeneous parts induced by ST, respec- 
tively. That is, we assume (p- m ii) = 0, where (■ • •) = 
J dr(- ■ -)/V represents the space average. Then p\ is re- 
lated to the space average of 6T by 



- (dp/dT) p (6T)(t), 



(2.4) 



which follows from the thermodynamic relation dp = 
(dp/dT) p dT + (dp/dp)rdp and the condition that the 
space average of the density deviation vanishes ((5p) = 
0). It is important that the combination p(r,t) + p c gz 
is nearly homogeneous or |pi(t)| ^S> |pinh(*"j t)\ for fluid 
motions much slower than the acoustic time t ac = L/c 
(c ~ 10 4 cm/s being the sound velocity) 24 

Now we derive the equation for ST from the heat con- 
duction equation 



PT 



d_ 

dt 



u-V )s = XV 2 ST, 



(2.5) 



where s(r,t) is the entropy per unit mass. Here Ss con- 
sists of the equilibrium part s cq (z) with 



and the nonequilibrium deviation 
Ss(r,t) = T~ X G P 



(2.6) 



(2.7) 



With the aid of the thermodynamic identity (dT/dp) s = 
(dT/dp) p (l~l/-f s ), we rewrite (2.5) to obtain the desired 
equation for ST, 



d_ 

dt 



v-V-DV 2 }ST= -a s v 



a s -(ST), (2. 



where D = X/C p is the thermal diffusivity and 

u s = 1 - 77 1 - 
The specific heat ratio 7 S behaves as 



C p /C\ 



-j+a 



>1, 



(2.9) 



(2.10) 



where C' p ~ e -7 and Cy ~ e~ a are the specific heats 
(per unit volume) at constant p and V, respectively, with 
a = 0.1. The first term on the right hand side of (2.8) 
arises from ds cq /dz. Inside plumes the temperature is 
adiabatically cooled if they go upward (v z > 0), or adi- 
abatically warmed if they go downward (v z < 0). In 
this way this term suppresses upward motion of warmer 
plumes from the bottom and downward motion of cooler 
plumes from the top, resulting in the Schwarzschild crite- 
rion of convection onset (the adiabatic temperature gra- 
dient effect). On the other hand, the second term arises 
from pi(t), leading to the piston effect (^4|. It is worth 
noting that the space integral of (2.8) in the cell becomes 



VC V — (ST) = A / dan- VST, 
dt J 



(2.11) 



where use has been made of (v) = 0. The right hand 
side represents the rate of heat supply from the bound- 
ary surface, where da is surface element and n is the 



outward surface normal. Its time-integration is the total 
heat supply expressed as Vp(Ss), resulting in 



C v (5T)(t) = p(6s)(t), 



(2.12) 



which also follows (2.4) and the space average of (2.7). 
The appearance of CV in on the left hand side of (2.12) is 
a natural consequence under the fixed volume condition. 
Notice that (2.7) can also be written as 



ST(r,t) = ^Ss(r,t)+pT 



1 

Cv 



1 

a, 



(6s)(t). (2.13) 



This relation holds even in gravity if Ss is the deviation 
of s — s eq (z) as in (2.7). In addition, the density deviation 
Sp = p — (p) is written in our approximation as 

Sp = pK T g(z - L/2) - pa p (ST - (ST)), (2.14) 

where Kt — (d p / dp)r / p and we have set (Sp) = 0. 

Since C p 3> Cy near the critical point, the homoge- 
neous part of ST (second term) in (2.13) can easily dom- 
inate over the inhomogeneous part (first term) even when 
Ss is localized near a heated wall. Indeed, if a thermal 
disturbance is produced within a thermal boundary layer 
with thickness t near the boundary, the ratio of the ho- 
mogeneous part (ex (Ss)) to the localized inhomogeneous 
part (oc Ss) in (2.13) is of order (7^ — 1)£/L where L is 
the characteristic system length. Temperature homog- 
enization is achieved when (j s — 1)£ 3> L. By setting 
£ = (Dti) 1 / 2 we obtain the time constant of this thermal 
equilibration (the piston time) in the form, 



h = L 2 /D( ls - If 



(2.15) 



Next we consider the momentum equation for the ve- 
locity field v(r,t). On long time scales, sound waves 
decay to zero and the incompressibility condition 



V v = 



(2.16) 



becomes nearly satisfied (^> ihomo) |@- Then the dissi- 
pation of v is produced by the shear viscosity 77 and the 
usual Navier-Stokes equation in the Boussinesq approxi- 
mation may be set up in the form |lj , 



9 ^ Pinh 

— + vV)v = -V 

dt J p 



a p gSTe z + -V 2 v, (2.17) 



where the inhomogeneous part p-^ ensures (2.16), e z is 
the unit vector along the z axis, and p is the average 
density. The two equations (2.8) and (2.17) are our fun- 
damental dynamic equations closed under (2.16). In the 
conventional theory jl]]|], (2.17) has been used, but the 
right hand side of (2.8) vanishes. 

As another characteristic feature near the critical 
point, the Prandtl number behaves as 



Pr = rj/pD ~ e 



(2.18) 



For example, Pr = 350 at T/T c - 1 = 10" 
This means that the time scale of the thermal diffusion 
is much slower than that of the velocity in the critical 
region. Based on this fact, the simulation in Ref. ]3l| ] 
was performed using the Stokes approximation in which 
the left hand side of (2.17) is set equal to zero. Good 
agreement with the experiments |p7| was then obtained 
for Ra co "/Ra c - 1 < 5 at e = 0.05. 

For Pr ^ 1, let us estimate the upper bound of Ra c °" 
below which the Reynolds number Re is smaller than 1 
or the Stokes approximation is allowable. The character- 
istic temperature variation (ST)± changing perpendicu- 
larly to the z axis and the characteristic velocity field v p \ 
are related by 



(2.19) 



where k ~ 2ir/L for roll patterns. If Ra co " / Ra c is con- 
siderably (but not much) larger than 1, (5T)±/AT is of 
order 1 (but somewhat smaller than 1). Then we obtain 

w„l ~ (Ra cmT /Ra c )D/L. (2.20) 

Thus the small Reynolds number regime is written as 

Ra co "/Ra c < Pr, (2.21) 

where use has been made of Re ~ v p \Lp/rj. For Pr 1 
there is a sizable range of Ra c °" in which the Stokes ap- 
proximation is justified. In passing, for < Ra COTr / Ra c — 
1 <C 1, the theory of the amplitude equation 6§] predicts 



v pl L/D ~ (ST) ± /AT ~ (Ra COII /Ra c - 1) 



1/2 



(2.22) 



from which we have Nu — 1 ~ Ra c °" / Ra c — 1 because 
the convective heat current is of order C p (ST)±v p i. In 
the next section we will estimate v p \ for much larger Ra. 

Analogously to (2.19), the inhomogeneous pressure de- 
viation pi n h is estimated as p; n h ~ (ot p p c g/k)(ST)x- If 
we assume Pl (t) ~ (dp/dT) s AT from (2.4) and AT ~ 
(ST)± as in (2.20), we find that p- m h/pi(t) is of order 
e~ 7 a g /T c fc and is much smaller than 1 from (2.2). This 
estimation justifies the assumption of the homogeneity of 
Sp(r,t) + p c gz made below (2.4). 



B. Free energy and heat production rate 

In the presence of small deviations of the temperature 
and the density, ST and Sp, around an reference equilib- 
rium state, we have an increase of the free energy func- 
tional SF. Up to the bilinear order of the deviations, it 
is of the form MM, 



SF = dr 



■(M 2 +9 z ^p 



(2.23) 



where the third term is the potential energy in gravity. 
All the deviations are assumed to change slowly in space 
compared with the thermal correlation length £. If we 
express Sp in terms of ST as in (2.14), we obtain 



SF 



1 

2T 



dr 



C P (ST - (ST)) 2 + C v {5Tf 



(2.24) 



where the constant term is omitted. We notice that SF 
decreases dramatically for j s ^> 1 in the process of adia- 
batic temperature homogenization. Furthermore, in the 
presence of velocity field, the total free energy change is 
the sum of SF and the kinetic energy of the velocity field, 



F K 



drpv z 



(2.25) 



Its time derivative is calculated from our dynamic equa- 
tions (2.7) and (2.17) in the form, 

f t (5F + F K ) = ~fdr(e th + e vis ) 

+ XT- 1 J da[ST(n-VST)}, [ ' 

where e t h and e V is are the thermal and viscous heat pro- 
duction rates (per unit volume) j34j ] , respectively, defined 
by 



6t h = XT-^WSTl 2 , 



rj^2(dvi/dxj) 

ij 



(2.27) 



(2.28) 



In the second term of (2.26) the surface integral is over 
the boundary of the cell, n being the outward unit vec- 
tor. In terms of the heat flux from the bottom Q, it is 
expressed as VQAT/TL if the top temperature is fixed. 



C. Basic relations in steady states 

We consider steady convective states in the Rayleigh- 
Benard geometry, in which the flow pattern is either time- 
independent not far above the convection onset or chaotic 
at larger Ra. We treat AT as a constant parameter. 
Under the condition of fixed heat flux at the bottom, 
however, AT(t) exhibits rapidly-varying fluctuations in 
chaotic states. In this case AT in the following relations 
represents the time-average of AT(t). The steady state 
averages (over space and time) will be denoted by (• • -) s 
to distinguish them from the space averages (• • •) used so 
far. 

We make (2.8) and (2.17) dimensionless by measur- 
ing space and time in units of L and L 2 /D and setting 
r = L~ 1 r and t = DL~ 2 t. The temperature deviation is 
written as 



ST(r,t)/AT= l-z + Ra,- 1 F(r,t), 



(2.29) 



where z = z/L. The dimensionless function T becomes 
nonvanishing in convective states and obeys 



Ra COII V z 



<*.-g{F) t (2.30) 



where V = L\7 is the space derivative in units of L. 
Then the (average) heat flux at the bottom is written as 
Q = (XAT/L)[1 + Ra-ifx], where 



A = -<(aF/a*)i=o>. 



The f\ is a function of Ra c °" and Pr. The Nusselt num- 
ber Nu — QL/XAT is expressed as 



Nu = 1 



Ra f\. 



(2.32) 



As the boundary condition of T we require T = at 
z = and 1 if T top and T^ ot are fixed. However, If T top 
and Q at the bottom are fixed, we have T = at z = 
and dT jdz = Ra{Nu — 1) at f = 0. The dimensionless 
velocity V(f,f) = (L/D)v obeys 

i- f JL + V • VJV = -VP inh + Te z + V 2 V, (2.33) 

where Pi„h ensures V • V = 0. 

Here we assume that the piston term, the second term 
on the right hand side of (2.30), can be neglected in 
steady states. For e = 0.05, the piston term in steady 
states is less than a few percents of the convection term 
v ■ V.F in (2.30) except at the boundaries. It thus pro- 
duces no significant effects on steady state heat transport 
(on Nu), while it can be crucial in the initial transient 
stage pl[ . Then, if the piston term in (2.30) is neglected, 
(2.30) and (2.33) become of the same form as those of 
usual incompressible fluids except that Ra COTT appears in 
place of Ra. At much smaller e, however, this assumption 
is questionable, because the noise part of (T) grows as 
e — * 0, as will be discussed later in the next section. We 
may conclude the following (at least for e = 0.05). (i) It 
follows the Gitterman-Steinberg criterion Ra COTI > Ra c 
in convective states in the compressible case p2] , ^3| . (ii) 
It is more nontrivial that the combination 



Ra(Nu-l) = fx(Ra co ", Pr) 



(2.34) 



should be a universal function of i?a corr and Pr from 
(2.32) in agreement with the experiments [ ^2|Jl^ ]. Notice 
that Ra(Nu — 1) = f\{Ra,Pr) holds for incompressible 
fluids in terms of the same fx- These experiments and 
more decisively that by Ahlers and Xu p5| indicate that 
fx should be nearly independent of Pr once Pr consider- 
ably exceeds 1. In the 3D simulation by Verzicco and Ca- 
mussi 1 41 , Nu became independent of Pr for Pr > 0.5. 



Theoretical support of this behavior using scaling argu- 
ments was presented in Ref. . 



In steady states we may also derive some exact rela- 
tions for variances among ST and v. Using the dynamic 
equations (2.8) and (2.17) we calculate the averages of 
d(6T) 2 /dt, dv 2 /dt, and d(z5T)/dt to obtain 



(|V£Zf 



+ a th (a th -a s )(Nu-l), (2.35) 



£ {{dv i /dx j f) s = Ra(D/L 2 ) 2 (Nu - 1). (2.36) 



(2.31) We a l so obtain a cross correlation, 



(v z 6T) s = a th D(Nu - 1), 



(2.37) 



which is nothing but the average convective heat flux (if 
XC P is multiplied). Here a th = AT/L = -(d6T/dz) s 
is the average temperature gradient and a g is the adia- 
batic temperature gradient defined by (1.2). If we use the 
usual hydrodynamic equations for incompressible fluids, 
the right hand side of (2.35) becomes a 2 h Nu, while (2.36) 
and (2.37) remain the same [§. In addition, (2.35) indi- 
cates ath > Qg in convective states in which Nu > 1. This 
is consistent with the convection criterion Ra c °" > Ra c . 
We obtain the averages of the two dissipation rates in 
(2.27) and (2.28) by multiplying X/T and 77 to (2.35) and 
(2.36), respectively. Using the thermodynamic identity 

ain 



Ta p — Cp(dT/dp) s , we obtai 



= T~ 1 Xa 2 h Nu, 



((eth) s " T 1 Aa t 2 h )/(e vis ) s = a th /a g 



(2.38) 



(2.39) 



The first relation (2.38) also follows from the average of 
(2.26). The second relation (2.39) holds only in convec- 
tive states (Nu > 1), while the right hand side is replaced 
by C p ath/Ta p pg = dth/% for the usual hydrodynamic 
equations of incompressible fluids. 



III. SIMULATION RESULTS 

We perform numerical analysis of (2.8) and (2.17) in 
2D using parameters of 3 He in a cell with L = 1.06 mm. 
The reduced temperature is e = 0.05 (except in Fig. 13), 
where j s = 22.8, Ta p = 26.9, A = 1.88 x 10" 4 erg/(cm 2 s 
K), D — 5.42 x 10~ 5 cm 2 /s, and Pr = 7.4 |l|,y,|l|. 
The condition (2.2) is well satisfied. The piston time t\ 
in (2.15) is given by 0.42 s. We apply a constant heat 
flux Q at the bottom z = for t > with a fixed top 
temperature T top at z = L. In steady states we have 
Ra cmT /Ra c = 0.90[AT/a g i - 1], where a g L = 3.57 uK. 
Thus (AT) on = 7.6 and Q on = 13.5 nW/s at the con- 
vection onset. We assume homogeneity of the boundary 
temperatures, T top and T^ ot , in the lateral x direction. 

In the experiments the aspect ratio was 57, so in the 
simulation |uj] the periodic boundary condition was im- 
posed in the x direction with period AL. This period 



was chosen because the roll period is close to 2L slightly 
above the onset for infinite lateral dimension |Q. Then, 
in steady states in the region 1 < Q/Q D n ^ 5, the linear 
relation 



Q/Qon - 1 = A [AT/(AT) oa - 1] 



(3.1) 



was numerically obtained with Aq = 2.2 in good 
agreement with the experiments. From Nu = 
[Q/AT]/[Qon/(Ar)on], the behavior of Nu is known 
from (3.2) in the range 1 < Q/Q on ^5 5. In particular, 
slightly above the onset, we have 



Nu — 1 = Ai(Ra co "/Ra c - 1) 



(3.2) 



where A\ = 0.64 in fair agreement with the theoretical 
value (Ai ^ 0.70 for Pr = 7.4) @. This behavior is also 
consistent with (2.22). 

In this work we are interested in fluid motion for rela- 
tively large Ra up to 3 x 10 6 . In the following we show 
two sets of the numerical results. In the first set, peri- 
odic sidewalls are assumed at x ~ and x — Lj_ with 
period L± = AL as in Ref. In Table 1 the steady 

state values of AT, Ra cmr , Ra, Nu, and Re are writ- 
ten, where Re is a Reynolds number to be defined in 
(3.11). They are obtained for Q = 0.0458 /iW/cm 2 (^ 
3.4Q on ), 0.965//W/cm 2 (= 71Q on ), and 122.2^W/cm 2 (= 
9 x 10 3 Q O n)- For the smallest Q the system tends to 
a time-independent convective state, as already studied 
in Ref. |n], while for the other values of Q the system 
tends to a chaotic state without macroscopic boundary 
shear flow. In the second set, we perform simulations 
for A = 1,2, and 3 with insulating and rigid sidewalls 
at x = and AL, at which v = and through which 
there is no heat flux (d5T/dx — 0), as will be presented 
in Figs.4, 12, and 13. 

In addition, if the temperature difference will be sim- 
ply written as AT, it should be taken as the time average 
of AT(t) in a steady state. We also assume that Pr is 
considerably larger than 1 in the following arguments. 



A. Transient behavior 

In Fig.l we show numerically calculated AT(t) = 
Tbot(i) - T p for Q = 0.965^W/cm 2 in (a) and for 
Q = 122.2/iW/cm 2 in (b). They nearly coincide with 
the upper broken curve without convection (v = 0) in 
the initial stage before the maximum is attained. The 
latter curve is calculated from (2.8) as 



[Ar«)]. = fV- 

A V 7T 



4- 



1 



(3.3) 



where t\ is defined by (2.15) and the integral in the brack- 
ets behaves as (irti/t) 1 / 2 for t 3> t\ E3]. If the piston 



term is absent and v — 0, (2.8) becomes the simple dif- 
fusion equation, yielding [AT(t)] = {2Q / ^(Dt/ir) 1 / 2 , 
which is about half of [AT(t)]o in (3.3) for t ^> t\ (see 
Fig. 3 in Ref. |H|]). We also show the numerically cal- 
culated AT(t) at fixed pressure where the piston term is 
absent (a s — in (2.8)) but v ^ 0. In (a) the experimen- 
tal curve is shown to have a lower peak and overdamp 
more slowly than in our simulation. In (b) the selected 
value of Q is in the region where no overshoot was ob- 
served in the experiment. See also Fig. 11, where the 
numerical curves of AT(t) will be given for other choices 
of the parameters. 

In Fig. 2 we show time evolution of the temperature 
profile at Q = 122.2 /iW/cm for periodic sidewalls. 
In (A) and (B) small-scale mushroom-like plumes are 
ejected from the bottom. In (C) and (D) they reach 
the top and are flattened there. In this initial stage the 
typical raising speed v p \ is estimated as L/t tI where t tT 
is the traversing time. From (A)-(C) we find that it is 
nearly equal to the free-fall velocity v g defined by 



= (Lga p AT) 1/2 = (RaPr) 1/2 D/L, 



(3.4) 



which is 2.37 cm/s. In this case the plumes leave the bot- 
tom at zero velocity and go upward with their velocity 
roughly of the form, 



1 - exp[-(t - t )/t vis ] 



(3.5) 



where to is the departure time, t v - ls ~ pR 2 jr\ is the vis- 
cous relaxation time with R being the plume size, and 



R 2 gpa p AT/r] 



(3.6) 



is the terminal velocity achieved by balance between the 
buoyancy and the viscous drag. For t tI <C i v is the viscous 
drag is negligible and we have v p i(t) ~ v 2 (t — t )/L and 
it r ~ L/vg. Thus, if the initial velocity is much less than 
v g , the free-fall condition becomes 



R/L > {Pr/Ra) 



1/4 



(3.7) 



under which Voo = (R/ L) 2 (Ra/ Pr) 1/2 v g > v g . In Fig.2, 
R/L - 1/3 and {Pr/Ra) 1 / 4 - 0.04, so the above condi- 
tion is satisfied. 

With the arrival of the plumes the heat current in- 
creases at the top, because T top is fixed, and a negative 
deviation of Ss is produced in a layer near the top. As can 
be known from (2.13), the piston effect is then operative, 
resulting in a homogeneous lowering of the temperature 
in the whole cell. In the time region around (E) the 
fluid is vigorously mixed with high Reynolds numbers. 
More precisely, the height-dependent Reynolds number 
Re(z,t) to be defined in (3.12) below is about 20 ex- 
cept in the vicinity of the boundaries. A downward flow 
of cooler fluid regions is then produced from the top. 
In the steady state (F), the temperature deviation be- 
comes considerably smaller than in the transient states, 



and the localized boundary shear flows are produced be- 
tween outgoing and incoming plumes with thickness £ v 
much smaller than L. 

The overshoot is more clearly illustrated in Fig. 3, 
which displays the average of ST(x, z, f) taken in the x 
direction, 



ST(z,t) = 



L± dx 

—6T(x,z,t), 

Li 



(3.8) 



for the points (A), (C), (E), and (F) in Fig. lb. As a char- 
acteristic feature, the temperature in the interior consists 
of global changes due to the piston effect and bumps 
due to localized plumes. In (E) the cooler layer becomes 
thicker temporarily near the top due to the excess heat 
flow. 

In our simulation the raising plumes leave the bottom 
and reach the top nearly simultaneously, resulting in a 
homogeneous temperature change, (i) Not far above the 
onset this mechanism is the main cause of the overshoot 
in compressible fluids. Note that a small peak appears in 
AT(t) even in the fixed pressure case (a s — 1) as shown 
in Fig. 2 of Ref. |3l) and as was observed by Behringer 
and Ahlers ^0|. Furthermore, in Ref. |3l[], the time scale 
of the overshoot (from the maximum to the minimum of 
AT(t)) due to the piston effect was predicted to be of 
order t D / (Ra co ™ / Ra c - 1), where t D = L 2 /AD{= 50 s) 
is the diffusion time. This fairly agrees with later anal- 
ysis of the experimental data [|TJ. (ii) For much larger 
Q such as those in Figs. la and lb, however, the down- 
ward flow from the top is also rapid enough to produce 
large overshoot, as demonstrated by the curves at fixed 
(height-dependent) pressure. Whether fixed is the vol- 
ume or the pressure, the time scale of the overshoot is of 
the order of the traversing time L/v s of the plumes due 
to gravity. 

As regards the overshoot behavior of AT(t), agreement 
between our simulation and the experiment |17|] becomes 
worse with increasing Q. We point out the possibility 
that in the experiment a synchronous arrival of plumes 
at the top might have not been realized for very large Q 
or for very short L/v g because of large lateral dimensions 
of the cell used. That is, if some plumes arrive at the top 
and others leave the bottom at the same time, negative 
interference between currents up and down will suppress 
overshoot. 



B. Steady state behavior 

Now we discuss the Nusselt number Nu in steady 
states. Fig. 4 shows the combination Ra COTI (Nu — 
l)/(Ra cmT ~ Ra c ) vs Ra co "/Ra c ~ 1 for periodic side- 
walls and for A = 1,2, and 3. This combination de- 
pends on Ra corr and A from (2.34) in steady states. The 
data (solid line) excellently agree with the numerical 



results for periodic sidewalls. We find that the scaling 
relation (1.5) nicely holds for Ra colT / Ra c > 10 for peri- 
odic sidewalls, while it holds only for Ra co " / Ra c > 10 3 
at A = 1. The exponent a in (1.5) is close to 2/7, but 
a = 1/4 is also consistent with our numerical data. If 
A ~ 1 and Ra c °" is not very large such that the plume 
size is of order L, large-scale fluid motions are suppressed 
by the rigid sidewalls. This marked tendency of the A 
dependent crossover of Nu was already reported in mea- 
surements for A = 0.5, 1, and 6.7 p0| |. It was also con- 
firmed in the 3D simulation by Kerr p2| for periodic 
sidewalls. 

In Fig. 5 we show the steady-state temperature devi- 
ation 8T(z) averaged in the x direction as in (3.8) and 
in time for the three values of Q in Table 1 for periodic 
sidewalls with period L± = 4L. The averages taken along 
the x direction become only weakly fluctuating in time 
in steady chaotic states (the relative fluctuations being 
of order 10% for the largest Q). As has been observed 
ubiquitously in the previous simulations, the tempera- 
ture gradient becomes localized within thermal bound- 
ary layers with thickness £t- Because AT = 2£tQ/X for 
It <C L, it is related to Nu by (1.6). The arrows in Fig. 5 
represent the maximum points, z — £ v and L — £ v , of the 
variance of the horizontal velocity defined by 



^-v x (x,z,ty 
L>± 



1/2 



(3.9) 



In Fig. 6 we plot the normalized velocity variances, 
v*(z)/v g in (a) and v*(z)/v g in (b), where v g is defined 
by (3.4) and 



dx 

l7 



v z (x,z.ty 



1/2 



(3.10) 



The time average of v 2 and v\ in the brackets is also 
taken in these figures. On one hand, v* take maxima at 
z = £ v and L — £ v , where £ v is hardly distinguishable 
from £t- On the other hand, v* is largest at the middle 
of the cell. We also find that the sum (the kinetic-energy 
variance) (v*) 2 + (v*) 2 is nearly constant in the interior, 
which was a finding reported in Ref. [^] . At large Ra the 
maxima of v* and v* are of the same order and will be 
identified as the typical plume velocity v p \. In our simu- 
lation we have v p i ~ 0.1u g (cx Ra 1 / 2 ), which is consistent 
with velocity measurements ||^2| . 

Kerr and Herring 0] made similar plots of the height- 
dependent velocity variances in their 3D simulations for 
free-slip sidewalls. They found that the characteristic 
length £ v defined by the peak positions of v*(z) becomes 
longer than £t — L/2Nu with increasing Ra; for exam- 
ple, for Pr = 7 they obtained £ v /£t ~ 1 at Ra — 10 4 and 
£v/£t ~ 3 at Ra — 10 7 . Verzicco and Camussi obtained 
a similar slow growing of £ v /£t at large Ra for Pr > 1 
in their 3D simulation with A = 1 [^l[. Also similarly, 
our 2D simulation with Pr — 7.4 gives £ v /£t = 2.54 and 



1.1 for Q = 122.2 and 0.965 /iW/cm, respectively, but we 
cannot draw a definite conclusion because of our limited 
range of Ra. 

In Fig. 7 we plot an overall Reynolds number Re vs 
R co " /R c — 1 in the simulation for periodic sidewalls. It 
is defined by 



Re 



(\v ■ Vv\ 



2„.|2\ 



1/2 



(3.11) 



where the averages are taken in the whole space region. 
The Re is smaller than 1 for R cmr /R c < 5 @. For 
larger values of R co ", it exceeds 1 and the effective ex- 
ponent <9(lni?e)/<9(ln i?a corr ) is from 1/4 to 1/3. How- 
ever, as suggested by Fig. 6, the strength of the veloc- 
ity fluctuations strongly depends on the distance from 
the boundary, so it is more informative to introduce a 
height-dependent Reynolds number, 



Re{z) 



dx\v ■ Vi>| 



dx\V 2 v\ 



-l 1/2 



(3.12) 



where the time averages of the integrands are taken. As 
shown in Fig. 8, Re{z) takes maxima at z ~ £ v and L — £ v 
of order 



Re(£ v ) ~ £ v v p ip/i], 



(3.13) 



where v p \ ~ v*(£ v ). This relation indicates Re(£ v ) ~ 
Ra 1 / 2 -' 1 with a ^ 2/7 from u pl - 0.1v g and £ v ~ £ T . The 
Re(z) becomes considerably smaller in the interior than 
at z ~ £ v , whose origin is the sparseness of the plumes 
in the interior (see (3.20) below). We confirm that Re 
is of the order of the space average dzRe(z)/ L. In 
the literature however, the (large-scale) Reynolds 

number has been identified as Re — v v \Lp/rj, which is 
much larger than Re(£ v ) in (3.13) by L/£ v . (For roll pat- 
terns, as was discussed below (2.21), we uniquely have 
Re = VpiLp/r).) 

At very large Ra the boundary layers should gradu- 
ally crossover from a laminar state to a turbulent state 
except within thin viscous sublayers with thickness Zq 
much shorter than £ v . In the inertial region zq < z < £ v 
of the boundary layer, it is natural to expect the loga- 
rithmic velocity profile p4|, 



v* x (z) = b^(a /p) 1/2 [Hz/z )+c l 



(3.14) 



where Co is the amplitude of the shear stress at the 
boundary with bo and cq being dimensionless numbers 
of order 1. We may set cto = r/lim^o D xz (z), where 
D X z{z) is the variance of the velocity gradient, 



D xz {z) = 



L dx ( d . . 
o L± \dz 



1/2 



(3.15) 



Then v* x (z) ! 
define zq by 



(ao/rj)z as z — > 0. It is appropriate to 



z o = v/(pvo) 1/2 , 



(3.16) 



which ensures Re(zo) ~ 1. The size of <7q should be equal 
to the typical size of pv x v z fit z — £y even if we consider 
localized shear flows, so we also have 



00 



(3.17) 



The ratio of the two lengths zq and £ v is given by 

£ v /z ~ v p i£ v p/r) ~ Re(£ v ), (3.18) 

which grows with increasing Ra. In Fig. 9a, v*(z) is fitted 
to the above logarithmic form in the inertial region for 
Q = 122.2 /iW/cm, where (<r /p) 1/2 = 0.067w g = 0.16 
cm/s, b = 1.2, c = 0.97, and z Q = 0.025L. In Fig.9b, 
we plot v*(z) and zD xz (z) on a logarithmic scale. We 
may conclude that these quantities do not behave as z in 
the inertial region of the boundary layers, although our 
Ra is not large enough to unambiguously demonstrate 
the logarithmic velocity profile. Here we point out that 
our results are not consistent with Shiraiman and Siggia's 
primary assumptions of £t < £ v and the linear profile of 
the mean shear flow, v x cx z, in the region z < £t 

In contrast to the averages taken along the x direc- 
tion, those taken along the z direction are rapidly varying 
functions of time at large Ra due to the random plume 
motions. We consider the vertical velocity variance de- 
fined by 



v*(x,t) 



dz 



v z {x,z,ty 



1/2 



(3.19) 



In Fig. 10 we display snapshots of v*(x, t), where the time 
average is not taken and peaks arising from the plumes 
become more apparent with increasing Q. For our Ra re- 
alized, the space regions occupied by the plumes become 
more sparse with increasing Ra in the interior. As the 
plumes move through the cell, they remain distinguish- 
able from the ambient fluid because the thermal diffusion 
length (DL/vpi) 1 ^ 2 does not much exceeds £ v . So we may 
define the volume fraction of the plumes <p p i. The convec- 
tive heat current is of order <fip\v p \CpAT ~ XNuAT/ L, 
leading to 



cj) p \ - D/£ T v ph 



(3.20) 



which is of order Pr^ 1 Re^y)^ 1 <C 1 from (3.13). For 
much larger Ra, the plumes will generate smaller scale 
eddies, ultimately leading to fully developed turbulence 
in the interior, as will be discussed in Section 4. 



C. Overall Temperature Fluctuations 

When a plume with a volume Vr reaches the bound- 
ary, it transfers a heat of order CpATVa to the boundary 
wall. As indicated by (2.13), the piston effect then gives 
rise to a homogeneous change in (ST) (t) of order 

(ST)n^i s {Va/V)AT. (3.21) 

Of course, the real plumes are extended objects and are 
continuously arriving at the boundary in high Ra con- 
vection. Thus Vpi/V in the above formula should be re- 
garded as the fluctuation amplitude of the plume volume 
fraction (j) p \ in the interior, although we do not know its 
dependence on Ra etc at present. If T top and Q at the 
bottom are fixed as in our simulation, AT(t) should also 
consist of fluctuations of the same origin. Because of the 
strong critical divergence of j s , we expect that the rela- 
tive amplitude (ST) a/ AT would increase as e is decreased 
with a fixed size of AT. 

Fig. 11 displays time sequences of (8T)(t) and AT(t) 
at fixed volume and pressure for periodic sidewalls with 
L±_ = 4T, which demonstrates strong correlations be- 
tween these two deviations at fixed volume. In case 
(a) (upper figure) we set e = 0.05 (7,. = 22.8), AT = 
0.17mK, Pr = 7.4, Ra = 7.38 x 10 4 , and Nu = 4.06, 
while in case (b) (lower figure) we set e = 0.01 (7 S = 119), 
AT = 0.19mK, Pr = 37.7, Ra = 4.14 x 10 5 , and 
Nu = 6.04. The steady state values of AT in the two 
cases are chosen to be only slightly different. At fixed vol- 
ume, the fluctuations of (5T)(t) and AT(t) are strongly 
correlated, and are larger and slower for (b) than for (a) 
in steady states (t > 100). This is because of the critical 
enhancement of the piston effect and the critical decrease 
of D with decreasing e. At fixed pressure, where the pis- 
ton effect is absent, AT(t) exhibits noises much smaller 
than those at fixed volume and (ST) (t) smoothly changes 
in time. It is worth noting that this noise increase at fixed 
volume accompanied with an increase of Ra is contrary to 
the usually measured noise behavior of the temperature. 
For non-critical fluids, if the temperature is measured at 
the center of a cell, its fluctuation amplitude divided by 
AT is known to decrease with increasing Ra as Ra~" n . 
The exponent /3„ was about 0.15 in a cell with A = 1 
ill- 

D. Random Reversal of Macroscopic Flow 

For a convection cell with A ~ 1, it is well-known 
that large-scale shear flow develops near the bound- 
ary of the cell for large enough Ra |l^,|2l],^2| . More- 
over, it has also been observed that the global circula- 
tion changes its orientation over long time scales plL[l9[ . 
For the case of A = 1, e = 0.05, Q = 40.7 /iW/ch?, 
Ra = 1.68 x 10 6 (^ Ra cml ), and Nu = 5.97, we plot a 



numerical time sequence of a circulation T(t) in Fig. 12. 
Here, 

F(t) = J f f dx[v x (x,L - d,t) - v x (x,d,t)]/L /^ 2 2) 
+ J d ~ dz[v z (L - d,z,t) - v z (d,z,t)]/L, 

where the integration is along a square contour with dis- 
tance d — 0.05L from the cell boundary. This quantity is 
positive for clockwise circulation and negative for coun- 
terclockwise circulation. In Fig. 12, AT(t) is also plotted, 
which exhibits particularly large fluctuations on the oc- 
casion of orientation changes. This is a natural result 
because large-scale reorganization of the flow pattern is 
needed for an orientation change. Fig. 13 illustrates the 
velocity patterns at t = 228,269, and 311s in Fig. 12. 
They closely resemble a picture of the measured velocity 
pattern in Ref. p2fl . 

IV. SCALING THEORY 

Rayleigh numbers realized in the existing simulations 
are still moderate in the sense that the plumes do not 
have enough kinetic energies such that they do not gen- 
erate fully developed turbulence in the interior. In this 
pre- asymptotic regime of steady states, we may under- 
stand the numerical and experimental data using a very 
simple zeroth- order theory. First, we set £ = It = £ v ne- 
glecting the possible small difference between £t and £ v 
mentioned below (3.10). The plume sizes in the horizon- 
tal direction are also of order £. Second, in our simulation 
the plumes are ejected into the interior with a velocity 
v p \, for which the viscous drag and the buoyancy are bal- 
anced or 

V r 2 v pl ~ ga p AT. (4.1) 

Thus v p \ is of the order of the terminal velocity Voa ~ 
RaD£ 2 /L 3 in (3.6) with R ~ £. In the interior we find 
that (i) gravity-induced acceleration of the plumes is sup- 
pressed by the viscous drag, (ii) (v*) 2 + (v* z ) 2 is nearly 
independent of z as stated below (3.10), and (iii) the last 
two terms on the left hand side of (2.17) are numerically 
of the same order. For example, the ratio of the average 
of (a p gST) 2 in the x direction to that of \(rj/ p)V 2 v\ 2 is 
about 4 at z ~ £ and is fluctuating around 1 in the inte- 
rior for the largest Q in Table 1. These support v p \ ~ 
in the interior. Third, to the sum rule (2.36) for the 
velocity gradients, the contribution from the boundary 
layers is of order v pl /£L, while that from the interior is of 
order <j> p \v pl /£ 2 ~ Dv p \/£ 3 from (3.20). If use is made of 
(4.1) and the sum rule (2.36), these boundary-layer and 
bulk contributions become both of order RaNu(D / L 2 ) 2 , 
which has also been confirmed numerically. Thus, 

v p i - Ra 1/2 D/L, (4.2) 



Nu ~ L/£ ~ l/0 p i~ Ra 1 ' 4 . (4.3) 

These quantities are independent of Pr. In particular, 
the independence of Nu on Pr is consistent with the ex- 
periments jjj|,[15],[l7j . Note that v p \ (~ Woo ) is smaller than 
v g in (3.4) by Pr- 1 / 2 . 

Our height-dependent Reynolds number at z = £ in 
(3.13) becomes 

Ee(€) ~ Ra 1/4 /Pr. (4.4) 

The usual large-scale Reynolds number is given by Re ~ 
VpiLp/rj ~ Ra 1 / 2 / Pr. As Re(£) exceeds a crossover value 
i?e*, plumes will induce turbulence in the interior. Our 
simple scaling theory is valid for Ra < (Re* Pr) 4 . In 
our simulation we have Re(£) = O.S&Ra 1 ' 4 , so that if we 
set Re* ~ 10 3 (regarding plumes as jets Q), the upper 
bound is crudely estimated as 5 x 10 13 . The transition 
from the scaling (4.3) to the asymptotic scaling occurs 
over a very wide range of Ra. Similarly, Grossmann and 
Lohse [|| considered a transition of a laminar boundary- 
layer flow to a turbulent boundary layer when the local 
Reynolds number on the scale of £ v at z ~ £ v exceeds a 
value of order 420. Then Nu was claimed to be better 
expressed by 

Nu ~ Ra 1/4 (l + dRa b ) (4.5) 

than the single power-law form, where C\ and b are small 
coefficient and exponent, respectively, dependent on Ra 
and Pr under investigation (both being of order 0.1). 
This proposed form of Nu was later claimed to be in 
good agreement with data (lj] . 

Here it would be informative to add more supplemen- 
tary explanations of the previous scaling theories, (i) Shi- 
raiman and Siggia assumed fully developed turbu- 
lence in the interior. Then the maximum of the turbulent 
velocity gradient is of order Sd = V^dl P — ( v piP/ Lv) 1 ' 2 
at the smallest eddy size kj 1 (~ (-q/ pvpiY^L 1 / 4 ) if the 
Kolmogorov cascade is assumed with the energy dissi- 
pation rate v^/ L [Q (the sparseness of ejected plumes 
being neglected). If the left hand side of the sum rule 
(2.36) is estimated as S 2 , it follows the relation, 

v p i ~ [PrNuRa) 1/z D/L. (4.6) 

Furthermore, they assumed the linear horizontal veloc- 
ity profile v x ~ (aa/rj)z in the region z < £t ~ L/Nu, 
where ao is given by (3.17). From the thermal diffusion 
equation v x dST/dx — D\7 2 ST (the time-dependent fluc- 
tuations being neglected), they obtained the scaling, 

£^ ~ fnfa/vDL, (4.7) 

by setting d/dx ~ L -1 and V 2 <~ (d/dz) 2 ~ i^ 2 for a 
cell with A ~ 1. From (4.6) and (4.7) they found 

Nu ~ Pr- 1/7 Ra 2/J . (4.8) 



However, as discussed below (3.14), our simulation sug- 
gests that the velocity deviates significantly from the lin- 
ear profile in the boundary layers, (ii) Castaing et al. || 
assumed the balance (4.1) at the length £t, 

v pl - £^ga p AT/T] - RaNu^ 2 D / L. (4.9) 

They furthermore assumed that the typical temperature 
scale in the interior is (8T) C ~ v^\ja v gL and that the av- 
erage heat current (= NuXAT/L) is of order C p (8T) c Vp\. 
From these relations {8T) C may be eliminated to give 
(4.6). If we combine (4.6) and (4.9), we are again led to 
(4.8). Therefore, to justify their arguments, the presence 
of fully developed turbulence in the interior seems to be 
required, (iii) Grossmann and Lohse estimated the 
bulk and boundary-layer contributions to the sum rules 
for the temperature gradient and the velocity gradients, 
the incompressible version of (2.35) and (2.36). Their 
primary assumption is that the boundary layer thickness 
for the velocity is given by £ v ~ L/Re 1 / 2 in terms of the 
large-scale Reynolds number Re p3. Note that this as- 
sumption is not consistent with our zeroth-order scaling 
theory with respect to the Pr dependence. In particular, 
in the case where Pr > 1 and the boundary-layer contri- 
butions are dominant both for the temperature and the 
velocity, they obtained Nu ~ Pr^ 1 / 12 Ra 1 / 4 . In this case 
we also find £ v /£t ~ Pr 1 / 3 from their theory. They pre- 
dicted that this pre-turbulent scaling crossovers to the 
asymptotic turbulent scaling very slowly as in (4.5). 



V. CONCLUDING REMARKS 

We have presented a hydrodynamic model of compress- 
ible fluids properly taking into account the piston effect 
and the adiabatic temperature gradient effect. Though 
performed in two dimensions, our simulation has revealed 
some new effects unique in near-critical fluids, such as 
the overshoot behavior and the amplification of the over- 
all temperature fluctuations as T — > T c . It generally 
explains the experimental findings Jl6|,[l^] , but a discrep- 
ancy remains in the overshoot behavior at high heat flux 
Q as discussed in Section 3. It is desirable to extend 
simulation to smaller e and higher Ra. Also more exper- 
iments on the overshoot and the temperature noises etc. 
are needed to resolve the discrepancy and to confirm the 
new predictions. As by-products, we have numerically 
examined steady state properties not treated in the pre- 
vious simulations, such as the logarithmic velocity pro- 
file and the random reversal of macroscopic shear flow. 
They are universal aspects present both in compressible 
and incompressible fluids. 

We have assumed that the fluid is in the supercritical 
region not very close to the critical point such that the 
conditions (2.1) and (2.2) are satisfied. However, if AT 



exceeds T — T c or if T top is below T c , we encounter a vari- 
ety of new effects such as boiling and wetting under heat 
flow and gravity (2j]j5^]. We believe that such problems 
should provide us a new challenging field in which nonlin- 
ear dynamics and phase transition dynamics are coupled. 
These problems are beyond the scope of this paper. 

We thank H. Meyer for valuable suggestions and com- 
ments. Thanks are also due to P. Tong for informative 
correspondence. This work is supported by Japan Space 
Forum grant H13-264. 



[1] S. Chandrasekhar, Hydrodynamic and Hydromagnetic 
Stability (Clarendon Press, Oxford, 1961). 

[2] E.D. Siggia, Annu. Rev. Fluid Mech. 26, 137 (1994). 

[3] B.I. Shiraiman and E.D. Siggia, Phys. Rev. A 42, 3650 
(1990). 

[4] LP. Kadanoff, Physics Today 54, No.8, 34 (2001). 

[5] B. Castaing, G. Gunaratne, F. Heslot, L.P. Kadanoff, A. 
Libchaber, S. Thomae, Xiao-Zhong Wu, S. Zaleski, and 
G. Zanetti, J. Fluid Mech. 204, 1 (1989). 

[6] S. Grossmann and D. Lohse, J. Fluid Mech. 407, 27 

(2000) . 

[7] S. Grossmann and D. Lohse, Phys. Rev. Lett. 86, 3316 

(2001) . 

[8] The usual hydrodynamic equations for one-component 
fluids in the Boussinesq approximation are valid only 
when (i) the specific heat ratio C p /Cv is close to 1, (ii) 
AT is much larger than a s L in (1.2), and (iii) the gravity- 
induced stratification is weak such that (2.2) holds. 

[9] M. Sano, Xiao-Zhong Wu, and A. Libchaber, Phys. Rev. 
A 40, 6421 (1989). 
[10] Xiao-Zhong Wu and A. Libchaber, Phys. Rev. A 45, 842 
(1992). 

[11] S. Ashkenazi, Ph.D. thesis, Weizmann Institute of Sci- 
ence, Rehovot, Israel, 1997(unpublished); S. Ashkenazi 
and V. Steinberg, Phys. Rev. Lett. 85, 3641 (1999). 

[12] X. Chavanne, Ph.D. thesis, Universitee Joseph Fourier, 
Grenoble, 1997 (unpublished). 

[13] X. Chavanne, F. Chilla, B. Castaing, B. Hebral, B. 
Chabaud and J. Chaussy, Phys. Rev. Lett. 79, 3648 
(1997); ; X. Chavanne, F. Chilla, B. Chabaud, B. Cas- 
taing, and B. Hebral, Phys. Fluids 13, 1300 (2001). 

[14] X. Xu, K.M.S. Bajaj and G. Ahlers, Phys. Rev. Lett. 84, 
4357 (2000). 

[15] G. Ahlers and X. Xu, Phys. Rev. Lett. 86, 3320 (2001). 
[16] A.B. Kogan, D. Murphy and H. Meyer, Phys. Rev. Lett. 

82, 4635 (1999). 
[17] A.B. Kogan and H. Meyer, Phys. Rev. E 63, 056310 

(2001). 

[18] J.J. Niemela, L. Skrbek, K.R. Sreenivasan and R.J. Don- 
nelly, Nature 404, 837 (2000). 

[19] J.J. Niemela, L. Skrbek, K.R. Sreenivasan and R.J. Don- 
nelly, J. Fluid Mech. 449 169 (2001). 

[20] B.J. Gluckman, H. Willaime, and J.PGollub, Phys. Flu- 
ids A 5, 647 (1993). 



[21] S. Cioni, S. Ciliberto, and J. Sommeria, J. Fluid Mech. 
335, 111 (1997). 

[22] X.-L. Qiu and P. Tong, Phys. Rev. E 64, 036304 (2001). 

[23] A. Onuki, Phase Transition Dynamics, (Cambridge Uni- 
versity Press, Cambridge, 2002). 

[24] A. Onuki and R.A. Ferrell, Physica A 164, 245 (1990); 

A. Onuki, H. Hao and R.A.Ferrell, Phys. Rev. A 41 2256 

(1990) . 

[25] J. Straub and K. Nitsche, Fluid Phase Equilibria, 88, 183 

(1993). J. Straub, L. Eicher and A. Haupt, Phys. Rev. E 

51, 5556 (1995). 
[26] H. Boukari, J.N. Shaumeyer, M.E. Briggs and RW. 

Gammon, Phys. Rev. A 41, 2260 (1990). 
[27] B. Zappoli, D. Bailly, Y. Garrabos, B. Le Neindre, and 

D. Beysens, Phys. Rev. A 41, 2264 (1990). 
[28] P. Guenoun, B. Khalil, D. Beysens, F. Kammoun, B. le 

Neindre, Y. Garrabos and B. Zappoli, Phys. Rev. E 47 

1531, (1993) 

[29] F. Zhong and H. Meyer, Phys. Rev. E 51, 3223 (1995) 
[30] S. Amiroudine, , P. Bontoux, P. Larroud, B. Gilly and 

B. Zappoli, J. Fluid Mech. 442, 119 (2001). 

[31] Y. Chiwata and A. Onuki, Phys. Rev. Lett. 87, 144301 
(2001). 

[32] M. Gitterman and V. Steinberg, J. Appl. Math. Mech. 
USSR 34, 305 (1971); M. Gitterman, Rev. Mod. Phys. 
50, 85 (1978). 

[33] P. Carles and B. Ugurtas, Physica D, 126, 69 (1999). 
[34] L.D. Landau and E.M. Lifshitz, Fluid Mechanics (Perg- 
amon, 1959). 

[35] E. E. DeLuca, J. Werne, and R. Rosner, Phys. Rev. Lett. 
64, 2370 (1990); J. Werne, E. E. DeLuca, R. Rosner, and 
F. Cattaneo, ibid. 67, 3519 (1991). 

[36] U. Hansen, DA. Yuen and S.E. Kroening, Phys. Flu- 
ids A 2, 2157 (1990); U. Hansen, D.A. Yuen and A.V. 
Malevsky, Phys. Rev. A 46, 4742 (1992). Here the Stokes 
approximation was used in the range 10 6 < Ra < 10 9 
in 2D. Nevertheless they could obtain the scaling (1.5) 
in agreement with the experiments, which suggests that 
the balance of the viscous drag and the buoyancy in (4.1) 
should have also been realized in the experiments. 

[37] S. Balachandar and L. Sirovich, Phys. Fluids A 3, 919 

(1991) . 

[38] S. Toh and E. Suzuki, Unstable and Turbulent Motion of 

Fluid, (World Scientific, 1993) p. 272. 
[39] C. Bizon, J. Werne, A. A. Predtechensky, K. Julien, W.D. 

McCormick, J.B. Swift and H.L. Swinney, Chaos 7, 107 

(1997). 

[40] S. L. Christie and JA. Domaradzki, Phys. Fluids A 5, 
412 (1993). 

[41] R. Verzicco and R. Camussi, J. Fluid Mech. 383, 55 

(1999) . 

[42] R. M. Kerr, J. Fluid Mech. 310, 139 (1996). 

[43] R. M. Kerr and J. R. Herring, J. Fluid Mech. 419, 325 

(2000) . 

[44] For < e 1 the scaling form a p — e~~<G(u) holds with 
u = (p/p c — For [w| < 1 or under (2.2) in gravity, 

we have G(u) G(0) and a v 2£ G(0)e" 7 . 

[45] G.P. Metcalfe and R.P. Behringer, J. Low Temp. Phys. 
78, 231 (1990). 

[46] To be precise, the time constant of the pressure- 
homogenization is given by thom = (iacii) 1,/2 under the 



condition t ac -C ii in terms of the acoustic time i ac = L/c 
and the piston time ti J23j . Here the damping of the pres- 
sure oscillation is mainly caused by the damping of the 
oscillatory heat current in the thermal diffusion layers. 
[47] L.D. Landau and E.M. Lifshitz, Statistical Physics (Perg- 
amon, New York, 1964), Chap. 12. It is well-known that 
the thermal fluctuations of 8T and 5p obey the Gaussian 
distribution proportional to exp(— 8F/kBT), where SF is 
given by (2.20). See Ref. [^3| for more discussions on this 
aspect. 



[48] G. Ahlers, M.C. Cross, P.C. Hohenberg, and S. Safran, 

J. Fluid Mech. 110, 297 (1981). 
[49] A. Schiiter, D. Lortz and F. Busse, J. Fluid Mech. 23, 

129 (1965). 

[50] R.P. Behringer and G. Ahlers, Phys. Lett. A 62, 329 
(1977). This experiment was performed in 4 He at T = 
2.184K at saturated vapor pressure. Here the fluid is 
nearly incompressible and 7 S = 1. 

[51] H. Meyer and A. B. Kogan, preprint. 

[52] A. Onuki, preprint. 



TABLE I. Parameters at e = 0.05 in steady states for periodic sidewalls 



Q (/xW/cm s) 


AT (mK) 


Ra co " 


Ra 


Nu-1 


Re 


0.0458 


0.0154 


3.43xlO a 


6.69 xW 


0.714 


0.655 


0.965 


0.135 


5.87 x 10 4 


5.54 x 10 4 


3.04 


3.035 


122.2 


6.89 


2.91xlO b 


2.91x10" 


9.29 


7.89 
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Fig.l. AT(t) vs time (solid line) calculated from (2.8) and (2.17) for (a) Q = 0.965 /iW/cm 2 and (b) Q = 122.2 
/iW/cm 2 . The temperature profiles for the points (□) on the curve in (b) are given in Fig. 2. The experimental data 
(+) are shown in (b). The upper broken curves in (a) and (b) represent the theoretical result (3.3) obtained 
from integration of (2.8) with v = 0. The dotted curves represent the numerical ones in the fixed pressure condition 
without the piston effect. 

Fig.2. Temperature profiles at (A), (B), (C), (D), (E), and (F) on the curve of Q = 122.2 ^W/cm 2 in Fig.lb (□). 
The temperature (and velocity) deviations are more enhanced in the transient states (A)~(E) than in a steady state 
(F). The ST at the bottom boundary z — is equal to AT(t) in Fig.lb. The plumes tend to be connected between 
bottom and top because Pr = 7.4. 

Fig.3. Time evolution of ST{z,t) defined by (3.5) at the points (A), (C), (E), and (F) in Fig.lb for Q = 122.2 
/^W/cm 2 . 

Fig. 4. Numerical results of Ra(Nu — l)/(Ra c °" — Ra c ) vs Ra corr /Ra c — 1 in steady states, obtained under the 
periodic boundary condition (+) and for A = 3(D), 2(*), and l(x). The first curve (+) is close to the experimental 
results for A = 57 @ (solid line) and is well fitted to the scaling form (1.5) with a ^ 2/7 for Ra co "/Ra c > 10 . With 
decreasing the aspect ratio A, crossover to the scaling occurs at much larger Ra COTr . 

Fig. 5. Height-dependent average temperature profiles 5T(z) divided by AT in steady states for the three Q values 
in Table 1. The arrows represent the maxima of v*(z) in Fig. 6a. 

Fig. 6. Normalized height-dependent variances, v*(z)/v s for the horizontal velocity in (a) and v*(z)/v s for the 
vertical velocity in (b) in steady states for the three Q values in Table 1. 

Fig. 7. Oveall Reynolds number Re defined by (3.11) as a function of Ra co "/Ra c — 1 in steady states for Q = 122.2 
fiW/cm 2 . 

Fig. 8. Height-dependent Reynolds number Re(z) defined by (3.9) in steady states for the three Q values in Table 

1. 

Fig. 9. (a) Height-dependent velocity variance v*(z) defined by (3.9) on a semi-logarithmic scale in steady states for 
Q = 122.2 fiW/cm 2 . (b) v*(z) (upperline) and velocity gradient variance zD xz (z) defined by (3.15) (lower line) on a 
logarithmic scale. 

Fig. 10. Snapshots of the normalized velocity variance v*(x, t)/vg averaged in the z direction defined by (3.17) for 
the three values of Q in Table 1. The system is periodic with period 4L in the x direction. The peak heights increase 
with increasing Q. For the largest Q this quantity changes in time as the plumes move in the cell, while for the other 
Q it is weakly dependent on or independent of time. 

Fig. 11. (5T)(t) and AT(t) at fixed volume (solid line) and at fixed pressure (broken line) for e = 0.05 (upper 
figure) and 0.01 (lower figure). The noises of these quantities at fixed volume increase as the reduced temperature e 
is decreased. 

Fig. 12. Time evolution of the circulation T(t) defined by (3.22) (upper figure) and AT(t) (lower figure) for Q = 122.2 
/iW/cm 2 in a cell with A = 1. The orientation of the macroscopic flow changes on a time scale of 50 s. The sign of 
F(i) represents the orientation of the macroscopic circulation, while the fluctuations of AT(t) become large when the 
orientation changes. 

Fig. 13. Velocity patterns at t = 228, 269, and 311 s for the run in Fig. 12. At t = 228 s the orientation is counter- 
clockwise, while at t = 311 s it is clockwise. At t = 269 s two large eddies with different orientations can be seen. 
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